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Abstract. - We modify the pre-factor of the semiclassical propagator to improve its efficiency 
in practical implementations. The new pre-factor represents the smooth portion of an orbit's 
contribution, and leads to fast convergence in numerical calculations. As an illustration of the 
accuracy and efficiency of the resultant propagator, we numerically calculate overlaps between 
quantum and semiclassical wave functions, as well as low-lying spectrum density in a 10- 
dimensional system contains unstable classical orbits. This sheds light on applying semiclassical 
propagator to high dimensional systems. 



Semiclassical propagator connects quantum dynamics with classical orbits. It enables one 
to evaluate quantum quantities from classical orbits. This is in principle applicable to high 
dimensional quantum systems as an alternative method to do first principle calculations. 
Since Gutzwiller adds Maslov phase to the Van Vleck semiclassical propagator [1], many tests 
confirm that the semiclassical propagator has remarkable accuracy [2-5]. Another appealing 
feature of the semiclassical method is easy to implement parallel computation by requiring 
each node to handle some orbits. The final result is simply a summation of each orbit's 
contribution. 

The fundamental approximation of semiclassical method is to treat the Hamiltonian locally 
as a quadratic function of coordinates and momenta. By expanding the Hamiltonian up to the 
second order around a phase space point, each orbit is associated to a time dependent quadratic 
quantum Hamiltonian in semiclassical calculations. In this sense, if numerical performance 
is not a concern, all kinds of semiclassical formulations are as accurate as the Van Vleck - 
Gutzwiller (VVG) propagator [1]. The later comes from the stationary phase approximation 
to the Feynman's path integral representation of the quantum propagator. In numerical 
calculations, however, the Herman and Kluk's (HK) formulation of semiclassical propagator [6] 
is often a favorite one [7-15]. 
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Like other formulations, the HK propagator suffers from the sign problem, i.e, fluctuations 
in nearby orbits' contributions. An orbit's contribution is usually a complex number with a 
phase related to the classical action of the corresponding orbit. The action changes rapidly 
from orbit to orbit. At the same time, the pre-factor of an orbit's contribution increases with 
time rapidly in the HK formulation, this makes the fluctuation increasing drastically with 
time, and it smears out any useful information sooner or later. The situation is especially 
worse for systems containing unstable orbits. 

A remedy to the sign problem is to use block method to smooth out the fluctuations [16-22] . 
The basic idea is same as the standard method to do numerical integration of a fluctuating 
function. Such treatment usually results in an exponential damping term that is equivalent 
to filter out contribution from unstable orbits. Another consequence of such treatment is 
that the formulation becomes complicated and the calculation of each orbit's contribution is 
numerically more costly. 

In this Letter, we present an approach to the sign problem by extracting smooth portion 
of an orbit's contribution. The resultant semiclassical propagator converges much faster while 
maintaining its original accuracy. We achieve this by modifying the pre-factor of the HK 
formulation with a positive factor. The semiclassical evolution operator with new pre-factor 
reads 

tsc{t) = / j2^W Rz ' tAz ^ SiZ ' t)/h |Z " A> (Z ' A| • (1) 

Here z = (q,p) stands for a 2 x TV dimensional phase space point with q and p being the 
coordinate and momentum respectively. An initial phase space point z arrives Zt — (qt,Pt) at 
time t along classical orbit. The action of the orbit is S(z,t) — J (pq — H)dt with H being 
the Hamiltonian of the system. The states \z, A) and \z t ,X) are Gaussian wave packets with 
A the squeeze parameter. In coordinate representation, it is 

/ ^ \ N/4, x ■ 

<a ^' A> = \^h) exp[ ~2h {x ~ q)2 + l h p{x ~ q)] - (2) 

The pre-factor i? Zjt is the square root of a determinant of the complex stability matrix, 

(3) 



R z ,t = \/det 



1 ( i 

- \M qq + M pp - i\M qp + -M pq 



where the TV x TV matrices M qq = d q qt(q,p), M pp = d p pt(q,p), M qp — d p qt(q,p), and M pq = 
d q pt(q,p) constitute the stability matrix M = d z zt{z) of the orbit. 

Without the modification factor A Z)t) Eq. Q is just the HK propagator. Note that the 
role of squeeze parameter A is equivalent to a scale transformation [19], i.e., one obtains the 
same result by applying the canonical transformation q — > V\q, p — > pj \f\ to the underlying 
classical dynamics while keeping A = 1 in to ©. In the following discussions, we set A = 1. 
With this setting, the modification factor has a simple form, 

A z>t = (2/3 + l) Ar [det(2/3M T Af + 1)]" 1/2 , (4) 

where (3 is a large positive number, and the superscript T stands for matrix transpose. Since 
the determinant of the symplectic stability matrix M is unit, det M = 1, and the inverse of M 
is already known, we use A a , t = (2/3 + l) Ar [det(2/3M T + M- 1 )]- 1 / 2 for practical calculations. 
Computation cost of the determinant of a 2N x 27V real positive matrix is about the same 
as that of a TV x TV complex matrix, thus calculation of the modification factor needs about 
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the same computer operations as that of the original pre-factor. This ensures the efficiency 
to implement the semiclassical evolution operator (JIJ. 

The modification factor J2J represents the smooth portion of an orbit's contribution. We 
obtain it from a kind of block treatment to the HK propagator. The starting point is similar 
to that in Refs. [16-20]. It is also referred as Filinov transformation to the HK propagator. 
The basic idea of this block approach is to replace one orbit's contribution with a weighted 
average of nearby orbits' contribution. For weight function in Gaussian form, the summation 
of the nearby orbits' contribution has analytic result. Here, the Hamiltonian within each 
block is approximated as a time dependent quadratic function of coordinates and momenta 
associated with a representing orbit. This is equivalent to treat the nearby orbits within a 
block via linearized dynamics associated with the representing orbit. It is apparent that the 
form of the weight function determines the performance of the resultant propagator. 

We perform the weighted average in coherent state (Gaussian wave packet) representation 
of HK propagator, (z' , X'\TJ^ C K (t)\zo, A'). Here XJ^ C K (t) is the semiclassical evolution operator 
l|T[l without the modification term A 2jt . In principle, any representation should give the same 
result, provided that one chooses proper weight function. The coherent state representation is 
over complete, it suffices to consider the case that the initial Gaussian is the same as the final 
one, Zq — zq. The squeeze parameter A' can be chosen arbitrary. We set A' = 1, i.e. equal to 
A in Eq. J3J). Under this specification, we choose a weight function in the following Gaussian 
form: 

W 9 (z - zT) = A w exp[- JrC* ~ z' t f + ie(z t - z' t )}. (5) 

Here Aw is the normalization constant. The phase points z and z' move to z t and z' t along 
classical orbits at time t, respectively. The width parameter /3 > is a large number to ensure 
that the weight function is effectively distributed near its center point. The imaginary part 
of the exponent is a small adjusted phase to cancel fluctuation in the average procedure. We 
choose e to make the orbit starting from z being a stationary phase orbit. 

After the weighted average, an orbit's contribution to the propagator gains an extra factor 

A z4 = ^ 2/3+ 1 ^ Ar y dz ' e -[2l3(6z t ) 2 + (Sz) 2 +2( Z - Z0 )Sz]/4h^ ^ 

where Sz = z' — z, Sz t = z' t ~ zt, and z is the initial phase space point of the orbit. We obtain 
the above result by applying linearized dynamics to the actions of the nearby orbits. Details 
of the procedure is similar to that in Refs. [16-20]. Another approximation is that the terms 
(Sz t ) 2 + 2(zt — Zo)5zt, as well as a small normalization term are dropped from the exponent 
of the integrand. These terms are small compared to the first term of the exponent. 

Note that the right hand side of JHJ is real and positive definite, i.e. the integrand's phase 
vanishes. This comes from two facts: (1) Since we set A' = A = 1, the second order terms of 
an nearby orbit's action cancel with the second order terms of the phase of the Gaussian wave 
packets; (2) The first order terms of the action plus the first order terms of the phase of the 
Gaussian wave packet are (pt — pa)Sqt — (qt — qo)Spt — (p — po)8q +{q — qo)Sp. Under linearized 
dynamics, (p - po)5q —(q— qo)Sp is symplectic, i.e., it is equal to (p t - p t)Sq t - (qt ~ qat)5p t , 
where zq — (qo,Po) moves to 2ot = (qotiPot) along classical orbit at time t. Thus the first 
order terms of the integrand's phase are linearly dependent on Szt = (Sqt,Spt). By proper 
choice of the adjust parameter e in the weight function (0, the integrand's phase vanishes. 

The integral in JBJ) represents an overlap of two classical density distributions at time 
t. One of them is initially p%(G) — exp{ — [(Sz) 2 + 2(z — Zo)Sz]/4h}. After evolution for a 
period of time t according to classical dynamics, its overlap with the density distribution 
P2 = exp[— f3(Sz) 2 /2h] gives the integral in JJJJ. 
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From consideration of classical dynamics, we argue that the second term in the exponent 
of pi is negligible. At initial time t = 0, the second density distribution p2 is highly localized 
in comparison with pi. Within the effective distributed area of P2, pi is virtually a constant. 
Thus p-2 dominates the integral and pi can be replaced by a constant at initial time. As time 
t increases, phase space points stretch out in some directions, and folding up in others. Such 
stretching and folding of phase space make the distribution pi distorted in phase space. The 
quadratic term of pi, i.e., pn = exp[— (5z) 2 /4fr] has only finite effective distribution area that 
is distorted during evolution. The distorted area overlaps only partially with p2 at time t. 
On the other hand, the second term of pi, i.e., P12 = exp[— (z — z )Sz/2h], however, behaves 
in a quite different way. This term is unbounded, and it is smoother than the quadratic 
term. Thus the distortion of phase space does not affect its overlap with the distribution p 2 . 
Similar to the initial time t — 0, within the effective distribution area of p 2 , P12 is effectively 
a constant, and can be dropped from the integral. 

As usual, we treat the overlap between pn and p2 at time t via linearized dynamics. The 
result is (@J. Since this modification factor is independent of the wave packet |zo, A'), we obtain 
the evolution operator in the form of Eq. Q . This concludes our theoretical derivation. 

We test the performance of the semiclassical evolution operator Q via Henon-Hcilcs (HH) 
Hamiltonian [23]. This system contains unstable classical orbits, and the original HK prop- 
agator converges only in a rather short time scale. Properties of the HH system is widely 
investigated classically, quantum mechanically, as well as semi-classically, including improve- 
ment to the semiclassical propagator by block treatments [17-19]. 

We first compare semiclassical wave function with exact quantum result via overlap be- 
tween semiclassical and quantum wave functions. Solid lines in Fig. ^ are the result of a 
particle of unit mass in the 2-dimensional Henon-Heiles potential v(x, y) = \{lo 2 x 2 + 0J 2 y 2 ) + 
Xy(x 2 + ijy 2 ), where ui x — 1.3, uj y — 0.7, A = —0.1, r\ = 0.1, and the Planck constant is 
set to 1. The wave function of the particle is initially a Gaussian with center position at 
(1, 1) and center momentum vanished. The semiclassical wave function is normalized to unit 
before doing overlap. These parameters are the same as that of Ref. [18], in which there is a 
similar calculation in a time scale between to 100, and it shows that the original HK prop- 
agator converges only for time less than 50. We see that, within a normalization constant, 
the semiclassical wave function possesses remarkable accuracy. As comparison, we show by 
dashed lines the corresponding result using formulation of Ref. [18], which is more suitable 
for calculation of wave function than other existing formulations. It is evident that for time 
t > 100, Eq. is more stable and converges faster. 

As is shown in Fig. ^ the accuracy of the semiclassical wave function is rather insensitive 
to the value of j3, provided that (3 is large enough, or the effective distribution area of the 
weight function is sufficiently localized. Small value of /3, or wider distributed weight function, 
leads to faster convergence, i.e., one needs less orbits. However, if (3 is improperly small, the 
linearized dynamics treatment to the modification factor is invalid beyond certain time. When 
this happens, we find that the normalization constant of the wave function becomes vanishing 
small after certain time. We use this property to check the validness of the width parameter 
p. 

Next, we test the semiclassical propagator's performance in finding energy spectrum of 
quantum systems. In fact, calculation of energy spectrum in high dimensional systems is 
fundamental important, and is still a challenge in many systems. Similar to Ref. [19], we 
make numerical test in a high dimensional Henon-Heiles system, i.e., a particle of unit mass 
moving in a A-dimensional HH potential, V(x) = \ Yli=i + ^ Z^Li^i^jM-i+^i+i)) wnere 
A = 0.1, and r\ = —0.3. The initial wave function is a Gaussian wave packet whose central 
position is qi — 2 and central momentum vanishes, Pi = 0. Classical orbit in such region is 
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more unstable than that of Fig. ^ One can obtain energy spectrum from auto-correlation 
function by Fourier transformation, or more advanced filter diagonalization [24] . In Fig. ffia) , 
we show the auto-correlation function versus time in a 2-dimensional case, and in Fig0 (b) 
we show the real part of its Fourier transformation, namely, the spectrum density. Here the 
thick solid line is the exact quantum result, and the dashed line is semiclassical result from 
Eq. (JJJ. The thin solid line is the semiclassical result using formulation of Ref. [19], an 
improved version of Walton's formulation [17]. We see that in a remarkable time scale, the 
semiclassical auto-correlation function is very close to the quantum result. Such time scale is 
enough to resolve the spectrum, and the semiclassical spectrum density is able to reproduce 
quantum result very well. This is especially the case for the peaks of the spectrum density, 
whose positions correspond to the energy spectrum. 

The behavior of the auto-correlation function is about the same as that of the wave function 
as shown in Fig. ^ Our test shows that the spectrum density converges faster than the wave 
function, and the peaks of the spectrum density are not sensitive to the width parameter (3. In 
calculations of Fig. El (a), we use 10 4 orbits by setting the width parameter to 1, 000. However 
the spectrum density of Fig. 0(b) is resultant of 2, 000 orbits with the width parameter of 100. 
In fact, auto-correlation function in shorter time is more important to determine the position 
of the peaks of the spectrum density, and shorter time wave function converges faster than 
that of longer time, and is more insensitive to the width parameter. This explains the fast 
convergence rate of the spectrum density. Note that, although Eq. gives more accurate 
correlation function in longer times than that of Ref. [19], the two semiclassical results in Fig. 
12 yield almost the same spectrum density. This property of the spectrum density sheds light 
on applying semiclassical propagator to find energy spectrum in high dimensional systems. 

Our numerical calculation for high dimensional HH system confirms the above properties 
in the lowlying energy. The lowlying spectrum density converges much faster than the higher 
ones, and the accuracy is quite insensitive to the width parameter /3. Figure|31is semiclassical 
spectrum density versus energy in a 10-dimensional HH system. To confirm the convergence, 
we make calculations for different width parameters j3. The thick solid, dashed, and dotted 
lines arc for width parameter (3 — 100, 200, and 50 respectively. From Fig. |3J it is evident 
that the spectrum density for energy less than 30 is well converged. Small j3 leads to fast 
convergence without sacrificing the accuracy of lowlying spectrum. We use 20,000, 60,000, 
and 10,000 orbits to produce the solid, dashed, and dotted lines, respectively. In fact, for 
(3 = 50, 2, 000 orbits is almost enough to obtain the lowlying spectrum, i.e., the peak positions 
of the spectrum density. Again, the thin solid line is the result using formulation of Ref. [19]. 
Although it gives almost the same spectrum, result from Eq. Q has better resolution, which 
demonstrates better accuracy in longer times. 

In summary, we have introduced a modification to the pre-factor of the semiclassical prop- 
agator in the HK formulation. It represents the smooth portion of an orbit's contribution to 
the semiclassical propagator. Unlike other formulations, this modification factor is just an 
overlap between two classical density distributions. Our treatment to this classical overlap 
yields a concise expression of the modification factor which is positive definite. The resul- 
tant semiclassical propagator converges much faster while maintaining its original accuracy. 
This improvement to the convergence is prominent in the calculation of lowlying spectrum 
density for high dimensional systems. This result sheds light on the efforts towards applying 
semiclassical propagator as a practical tool of first principle computation. 
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Fig. 1 - Overlap between normalized semiclassical wave function and exact quantum wave function 
versus time. Solid and dashed lines are results of Eq. Q and Herman's formulation [18], respectively. 
is the width parameter, and n is the number of orbits used in semiclassical calculations. 
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energy 

Fig. 2 - (a) Real part of the auto-correlation function C(t) versus time, (b) Real part of the 
spectrum density versus energy. Thick solid, dashed and thin solid lines are exact quantum result 
and semiclassical result from Eq. Q, as well as semiclassical result using formulation of Ref. [19], 
respectively. 
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Fig. 3 - Spectrum density versus energy in 10-dimensional HH system. The number of orbits n and 
the width parameter f3 in semiclassical calculations are indicated accordingly. The thin solid line is 
the semiclassical result of Ref. [19] with /3 = 100 and n — 10000. 



